#install.packages("xlsx", "plyr", "data.table", "ggplot2", "RColorBrewer", "stringr", "ggpubr")

library(xlsx)
library(plyr)
library(data.table)
library(ggplot2)
library(RColorBrewer)
library(stringr)
library(ggpubr)

setwd("H://OneDrive//研究//1-全球变化本土外来竞争Meta//2022//数据与分析")  #设置数据读取与输出文件夹

##################
###  数据输入  ###
##################
rm(list=ls())
Data <- read.xlsx("20220523.xlsx",sheetName = "data", encoding = "UTF-8")  #读取待分析数据
M_GCFs <- data.frame(m = c("DC", "DN", "NC", "TC", "TD", "TN", "PN"),
                     s1 = c("D", "D", "N", "T", "T", "T", "P"),
                     s2 = c("C", "N", "C", "C", "D", "N", "N"))  #两GCFs组合关系

###图1（交互关系）数据，自定义
D1 <- list(
           data.frame(x = factor(c("CK", "A", "B", "A+B Additive", "A+B Synergistic", "A+B Antagonistic"), levels = c("CK", "A", "B", "A+B Additive", "A+B Synergistic", "A+B Antagonistic")),
                      y = c(20, 25, 32, 37, 41, 16)),
           data.frame(x = factor(c("CK", "A", "B", "A+B Additive", "A+B Synergistic", "A+B Antagonistic"), levels = c("CK", "A", "B", "A+B Additive", "A+B Synergistic", "A+B Antagonistic")),
                      y = c(20, 15, 8, 3, 1, 12)),
           data.frame(x = factor(c("CK", "A", "B", "A+B Additive", "A+B Synergistic", "A+B Antagonistic"), levels = c("CK", "A", "B", "A+B Additive", "A+B Synergistic", "A+B Antagonistic")),
                      y = c(20, 25, 8, 13, 16, 7)),
           data.frame(x = factor(c("CK", "A", "B", "A+B Additive", "A+B +Synergistic", "A+B -Synergistic", "A+B +Antagonistic", "A+B -Antagonistic"), levels = c("CK", "A", "B", "A+B Additive", "A+B +Synergistic", "A+B -Synergistic", "A+B +Antagonistic", "A+B -Antagonistic")),
                      y = c(20, 25, 8, 13, 29, 3, 10, 18)),
           data.frame(x = factor(c("CK", "A", "B", "A+B Additive", "A+B Synergistic", "A+B Antagonistic"), levels = c("CK", "A", "B", "A+B Additive", "A+B Synergistic", "A+B Antagonistic")),
                      y = c(0, 5, 7, 12, 16, 4)),
           data.frame(x = factor(c("CK", "A", "B", "A+B Additive", "A+B Synergistic", "A+B Antagonistic"), levels = c("CK", "A", "B", "A+B Additive", "A+B Synergistic", "A+B Antagonistic")),
                      y = c(0, -5, -7, -12, -17, -6)),
           data.frame(x = factor(c("CK", "A", "B", "A+B Additive", "A+B +Synergistic", "A+B -Synergistic", "A+B +Antagonistic", "A+B -Antagonistic"), levels = c("CK", "A", "B", "A+B Additive", "A+B +Synergistic", "A+B -Synergistic", "A+B +Antagonistic", "A+B -Antagonistic")),
                      y = c(0, 5, -12, -7, 11, -17, -9.5, -1)))

###图3文献与研究，自定义
D3 <- list(
           read.xlsx("E:\\R\\working_area\\data.xlsx",sheetName = "literature",colIndex = c(1:5)), #文献国别 Literature country & 研究国别 Research country
           read.xlsx("E:\\R\\working_area\\data.xlsx",sheetName = "literature",colIndex = c(6:8)),  #文献、研究年份 Literature & Research year
           read.xlsx("E:\\R\\working_area\\data.xlsx",sheetName = "literature",colIndex = c(12:16)), #年份、国别 Literature & Research year & country
           read.xlsx("E:\\R\\working_area\\data.xlsx",sheetName = "literature",colIndex = c(12:14))  #文献、研究类型 Literature & Research type
)


############################
###  创建自定义统计函数  ###
############################

bootstrap <- function(data, tv, cv, group_c = "Target_Type", group_t, n_iter, n_samp) {  
                                               #data：需数据框格式, tv：处理变量, cv：对照变量, group_c：对照分组（默认物种）, group_t：处理分组, n_iter：迭代次数, n_samp：抽样数
  library(plyr)                                #dlply拆分数据框
  library(data.table)                          #data.table包tstrsplit、rbindlist函数
  
  group <- c(group_t, group_c)                 #分组变量
  
  data_c <- dlply(data, group_c)               #按组将数据拆分进list
  C_names <- names(data_c)                     #将分组组合
  C_g <- as.data.frame(tstrsplit(C_names, ".", fixed=TRUE))
                                               #分组组合拆为数据框
  names(C_g) <- group_c                        #命名分组数据框列名
  Mc <- list()                                 #对照均值list
  for (i in c(1:length(data_c))) {             #按分组数循环
    set.seed(1234)                             #设置抽样种子，防止多次运算结果存在差异
    S <- replicate(n_iter, sample(data_c[[i]][cv][,1], n_samp, replace = TRUE)) 
                                               #n_iter次迭代，放回抽样n_samp次的结果，矩阵S
    Mc[[i]] <- cbind(data.frame(C_g[i,]) ,apply(S, 2, mean, na.rm=TRUE), row.names = NULL)
                                               #对矩阵S求平均数，并添加分组变量作为list中的数据框
    names(Mc[[i]]) <- c(group_c, "v")          #给数据框命名
  }
  Mc <- as.data.frame(rbindlist(Mc))           #对照均值list转数据框
  
  data_t <- dlply(data, group)                 #按组将数据拆分进list
  T_names <- names(data_t)                     #将分组组合
  T_g <- as.data.frame(tstrsplit(T_names, ".", fixed=TRUE))
                                               #分组组合拆为数据框
  names(T_g) <- group                          #命名分组数据框列名
  Mt <- list()                                 #处理均值list
  for (i in c(1:length(data_t))) {             #按分组数循环
    set.seed(1234)                             #设置抽样种子，防止多次运算结果存在差异
    M <- replicate(n_iter, sample(data_t[[i]][tv][,1], n_samp, replace = TRUE))
                                               #n_iter次迭代，放回抽样n_samp次的结果，矩阵S
    Mt[[i]] <- cbind(data.frame(T_g[i,]) ,apply(M, 2, mean, na.rm=TRUE), row.names = NULL)
                                               #对矩阵S求平均数，并添加分组变量作为list中的数据框
    names(Mt[[i]]) <- c(group, "v")            #给数据框命名
  }
  
  E <- list()                                  #效应list
  for (i in c(1:length(Mt))) {                 #按分组数循环
    E[[i]] <- cbind(data.frame(T_g[i,]), Mt[[i]]$v - Mc[Mc[group_c] == as.character(unique(Mt[[i]][group_c])), "v"], row.names = NULL)
    names(E[[i]]) <- c(group, "E")             #给数据框命名
  }
  return(E)
}

mq <- function(data, v, group) {               #将data数据框按group分组计算v的平均值mean，2.5%和97.5%的分位数down、up，和大于0的比例p
  library(plyr)                                #dlply拆分数据框
  datac <- ddply(data, group,
                 .fun = function(xx, col) {
                   c(mean = mean     (xx[[col]], na.rm=T),
                     down = quantile (xx[[col]], 0.025, na.rm=T),
                     up   = quantile (xx[[col]], 0.975, na.rm=T),
                     p    = length(which(xx[[col]]>0))/length(xx[[col]])
                   )
                 },
                 v
  )
  names(datac) <- c(group, "mean", "down", "up", "p")
  datac <- datac[which(is.nan(datac$mean) == F),]
  return(datac)
}

replace <- function(vector, from, to) {       #批量替换向量vector中的元素，将from元素替换为to元素
  x <- vector
  for (i in c(1:length(from))) {
    x[x == from[i]] <- to[i]
  }
  return(x)
}

inter <- function(list, ab, a, b, v) {      #计算交互效应，list列表，ab、a、b交互、因子1、因子2序号向量，v比较的数据
  out <- list()
  for (i in c(1:length(ab))) {
    up <- max(quantile(list[[a[i]]][v], 0.975, na.rm=T), quantile(list[[b[i]]][v], 0.975, na.rm=T))
    down <- min(quantile(list[[a[i]]][v], 0.025, na.rm=T), quantile(list[[b[i]]][v], 0.025, na.rm=T))
    E <- as.vector(unlist(list[[a[i]]][v] + list[[b[i]]][v]))
    Emean <- mean(E, na.rm=T)
    Edown <- quantile (E, 0.025, na.rm=T)
    Eup   <- quantile (E, 0.975, na.rm=T)
  list[[ab[i]]]$inter[unlist(list[[ab[i]]][v]) - E > Edown - Emean & unlist(list[[ab[i]]][v]) - E < Eup - Emean] <- "Additive"
  if (mean(unlist(list[[a[i]]][v]), na.rm=T) / mean(unlist(list[[b[i]]][v]), na.rm=T) < 0) {
    list[[ab[i]]]$inter[((unlist(list[[ab[i]]][v]) - E) - (Eup - Emean)) * Emean > 0 & 
                    (unlist(list[[ab[i]]][v]) - up) * (unlist(list[[ab[i]]][v]) - down) > 0] <- "+Synergistic"
    list[[ab[i]]]$inter[((unlist(list[[ab[i]]][v]) - E) - (Eup - Emean)) * Emean > 0 & 
                    (unlist(list[[ab[i]]][v]) - up) * (unlist(list[[ab[i]]][v]) - down) < 0] <- "-Antagonistic"
    list[[ab[i]]]$inter[((unlist(list[[ab[i]]][v]) - E) - (Eup - Emean)) * Emean < 0 & 
                    (unlist(list[[ab[i]]][v]) - up) * (unlist(list[[ab[i]]][v]) - down) < 0] <- "+Antagonistic"
    list[[ab[i]]]$inter[((unlist(list[[ab[i]]][v]) - E) - (Eup - Emean)) * Emean < 0 & 
                    (unlist(list[[ab[i]]][v]) - up) * (unlist(list[[ab[i]]][v]) - down) > 0] <- "-Synergistic"
  } else {
    list[[ab[i]]]$inter[list[[ab[i]]][v] - E < Edown - Emean] <- "Antagonistic"
    list[[ab[i]]]$inter[list[[ab[i]]][v] - E > Eup - Emean] <- "Synergistic"
  }
  out[[i]] <- list[[ab[i]]]
  }
  return(out)
}

si <- function(data, v, group) {               #将data数据框按group分组计算v的平均值mean，2.5%和97.5%的分位数down、up，和大于0的比例p
  library(plyr)                                #dlply拆分数据框
  datac <- ddply(data, group,
                 .fun = function(xx, col) {
                   c(mean = mean     (xx[[col]], na.rm=T),
                     down = quantile (xx[[col]], 0.025, na.rm=T),
                     up   = quantile (xx[[col]], 0.975, na.rm=T),
                     p    = length(which(xx[[col]]>0))/length(xx[[col]])
                   )
                 },
                 v
  )
  names(datac) <- c(group, "mean", "down", "up", "p")
  datac <- datac[which(is.nan(datac$mean) == F),]
  return(datac)
}

replace <- function(vector, from, to) {       #批量替换向量vector中的元素，将from元素替换为to元素
  x <- vector
  for (i in c(1:length(from))) {
    x[x == from[i]] <- to[i]
  }
  return(x)
}


##################
###  统计分析  ###
##################

#效应值计算：E-Effect size，P-Performance，R-RCI
##全球变化因子
E_P_GCF <- bootstrap(Data, "ReTG", "ReTC", group_c = "Target_Type", c("GCFs", "NUM_GCFs_Record"), 1000, 100)
E_C_GCF <- bootstrap(Data, "TGRCI", "TCRCI", group_c = "Target_Type", c("GCFs", "NUM_GCFs_Record"),  1000, 100)
##全球变化因子数
#E_P_F <- bootstrap(Data, "ReTG", "ReTC", group_c = "Target_Type", "NUM_GCFs_Record", 1000, 100) 
#E_C_F <- bootstrap(Data, "TGRCI", "TCRCI", group_c = "Target_Type", "NUM_GCFs_Record", 1000, 100)
##固氮能力
E_P_N <- bootstrap(Data, "ReTG", "ReTC", group_c = "Target_Type", "NFIX_Target", 1000, 100) 
E_C_N <- bootstrap(Data, "TGRCI", "TCRCI", group_c = "Target_Type", "NFIX_Target", 1000, 100)
##生活周期
E_P_LC <- bootstrap(Data, "ReTG", "ReTC", group_c = "Target_Type", "Lifecycle_Target", 1000, 100) 
E_C_LC <- bootstrap(Data, "TGRCI", "TCRCI", group_c = "Target_Type", "Lifecycle_Target", 1000, 100)
##功能群
E_P_FG <- bootstrap(Data, "ReTG", "ReTC", group_c = "Target_Type", "FUN_Group_Target", 1000, 100) 
E_C_FG <- bootstrap(Data, "TGRCI", "TCRCI", group_c = "Target_Type", "FUN_Group_Target", 1000, 100)
##测定指标类型
E_P_I <- bootstrap(Data, "ReTG", "ReTC", group_c = "Target_Type", "Type_IND", 1000, 100) 
E_C_I <- bootstrap(Data, "TGRCI", "TCRCI", group_c = "Target_Type", "Type_IND", 1000, 100)

#计算差异值Difference size：仅GCF-M与GCF-S之间
Index <- unique(rbindlist(E_P_GCF)[,c(1,3)])
Index$Num <- c(1:nrow(Index))
Index_M <- Index[nchar(Index$GCFs) == 2,]
for (i in c(1:nrow(Index_M))) {
  Index_M$GCF1[i] <- substring(Index_M[i,1], 1, 1)
  Index_M$GCF2[i] <- substring(Index_M[i,1], 2, 2)
}
Index_M$Num1 <- as.numeric(replace(Index_M$GCF1, c("C", "D", "N", "P", "T"), c(1, 3, 9, 13, 17)))
Index_M$Num2 <- as.numeric(replace(Index_M$GCF2, c("C", "D", "N", "P", "T"), c(1, 3, 9, 13, 17)))
D_P_GCF <- list()
for (i in c(1:nrow(Index_M))) {
  D_P_GCF[[i]] <- cbind(E_P_GCF[[Index_M$Num[i]]][,c(1:3)],
                        data.frame(D = E_P_GCF[[Index_M$Num[i]]]$E - (E_P_GCF[[Index_M$Num1[i]]]$E + E_P_GCF[[Index_M$Num2[i]]]$E)))
  }
D_C_GCF <- list()
for (i in c(1:nrow(Index_M))) {
  D_C_GCF[[i]] <- cbind(E_C_GCF[[Index_M$Num[i]]][,c(1:3)],
                        data.frame(D = E_C_GCF[[Index_M$Num[i]]]$E - (E_C_GCF[[Index_M$Num1[i]]]$E + E_C_GCF[[Index_M$Num2[i]]]$E)))
}

#交互判定
I_P_GCF <- inter(E_P_GCF, 
                 c(5, 6, 7, 8, 11, 12, 15, 16, 19, 20, 21, 22, 23, 24), 
                 c(3, 4, 3, 4, 9, 10, 13, 14, 17, 18, 17, 18, 17, 18), 
                 c(1, 2, 9, 10, 1, 2, 9, 10, 1, 2, 3, 4, 9, 10), "E")
I_C_GCF <- inter(E_C_GCF, 
                 c(5, 6, 7, 8, 11, 12, 15, 16, 19, 20, 21, 22, 23, 24), 
                 c(3, 4, 3, 4, 9, 10, 13, 14, 17, 18, 17, 18, 17, 18), 
                 c(1, 2, 9, 10, 1, 2, 9, 10, 1, 2, 3, 4, 9, 10), "E")

#合并效应值、差异值
SE_P_GCF <- mq(rbindlist(E_P_GCF), "E", c("Target_Type", "GCFs", "NUM_GCFs_Record"))
SE_C_GCF <- mq(rbindlist(E_C_GCF), "E", c("Target_Type", "GCFs", "NUM_GCFs_Record"))
SD_P_GCF <- mq(rbindlist(D_P_GCF), "D", c("Target_Type", "GCFs"))
SD_C_GCF <- mq(rbindlist(D_C_GCF), "D", c("Target_Type", "GCFs"))
SI_P_GCF <- as.data.frame(table(rbindlist(I_P_GCF)$GCFs, rbindlist(I_P_GCF)$Target_Type, rbindlist(I_P_GCF)$inter))
names(SI_P_GCF) <- c("GCFs", "Target_Type", "inter", "Freq")
SI_C_GCF <- as.data.frame(table(rbindlist(I_C_GCF)$GCFs, rbindlist(I_C_GCF)$Target_Type, rbindlist(I_C_GCF)$inter))
names(SI_C_GCF) <- c("GCFs", "Target_Type", "inter", "Freq")
SE_P_F <- mq(rbindlist(E_P_GCF), "E", c("Target_Type", "NUM_GCFs_Record"))
SE_C_F <- mq(rbindlist(E_C_GCF), "E", c("Target_Type", "NUM_GCFs_Record"))
SD_P_F <- mq(rbindlist(D_P_GCF), "D", c("Target_Type", "NUM_GCFs_Record"))
SD_C_F <- mq(rbindlist(D_C_GCF), "D", c("Target_Type", "NUM_GCFs_Record"))
SE_P_N <- mq(rbindlist(E_P_N), "E", c("Target_Type", "NFIX_Target"))
SE_C_N <- mq(rbindlist(E_C_N), "E", c("Target_Type", "NFIX_Target"))
SE_P_LC <- mq(rbindlist(E_P_LC), "E", c("Target_Type", "Lifecycle_Target"))
SE_C_LC <- mq(rbindlist(E_C_LC), "E", c("Target_Type", "Lifecycle_Target"))
SE_P_FG <- mq(rbindlist(E_P_FG), "E", c("Target_Type", "FUN_Group_Target"))
SE_C_FG <- mq(rbindlist(E_C_FG), "E", c("Target_Type", "FUN_Group_Target"))
SE_P_I <- mq(rbindlist(E_P_I), "E", c("Target_Type", "Type_IND"))
SE_C_I <- mq(rbindlist(E_C_I), "E", c("Target_Type", "Type_IND"))

#获取样本量
N_P_GCF <- ddply(Data, c("Target_Type", "GCFs"),
                 .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_P_GCF)[3] <- "N"
N_C_GCF <- ddply(Data, c("Target_Type", "GCFs"),
                 .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_C_GCF)[3] <- "N"
N_P_F <- ddply(Data, c("Target_Type", "NUM_GCFs_Record"),
               .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_P_F)[3] <- "N"
N_C_F <- ddply(Data, c("Target_Type", "NUM_GCFs_Record"),
               .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_C_F)[3] <- "N"
N_P_N <- ddply(Data, c("Target_Type", "NFIX_Target"),
               .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_P_N)[3] <- "N"
N_C_N <- ddply(Data, c("Target_Type", "NFIX_Target"),
               .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_C_N)[3] <- "N"
N_P_LC <- ddply(Data, c("Target_Type", "Lifecycle_Target"),
                .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_P_LC)[3] <- "N"
N_C_LC <- ddply(Data, c("Target_Type", "Lifecycle_Target"),
                .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_C_LC)[3] <- "N"
N_P_FG <- ddply(Data, c("Target_Type", "FUN_Group_Target"),
                .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_P_FG)[3] <- "N"
N_C_FG <- ddply(Data, c("Target_Type", "FUN_Group_Target"),
                .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_C_FG)[3] <- "N"
N_P_I <- ddply(Data, c("Target_Type", "Type_IND"),
               .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_P_I)[3] <- "N"
N_C_I <- ddply(Data, c("Target_Type", "Type_IND"),
               .fun = function(xx, col) {length(xx[[col]])},"ReTG")
names(N_C_I)[3] <- "N"

#数据整理
P_F <- merge(N_P_F, SE_P_F, by = c("Target_Type", "NUM_GCFs_Record"), all = TRUE) 
C_F <- merge(N_C_F, SE_C_F, by = c("Target_Type", "NUM_GCFs_Record"), all = TRUE) 
P_F$GCFs[P_F$NUM_GCFs_Record == 1] <- "GCF(S)"
P_F$GCFs[P_F$NUM_GCFs_Record == 2] <- "GCF(M)"
C_F$GCFs[C_F$NUM_GCFs_Record == 1] <- "GCF(S)"
C_F$GCFs[C_F$NUM_GCFs_Record == 2] <- "GCF(M)"
P_GCF <- merge(N_P_GCF, SE_P_GCF, by = c("Target_Type", "GCFs"), all = TRUE) 
C_GCF <- merge(N_C_GCF, SE_C_GCF, by = c("Target_Type", "GCFs"), all = TRUE) 
P_GCF <- rbind(merge(N_P_GCF, SE_P_GCF, by = c("Target_Type", "GCFs"), all = TRUE), P_F)
C_GCF <- rbind(merge(N_C_GCF, SE_C_GCF, by = c("Target_Type", "GCFs"), all = TRUE), C_F)
names(SD_P_F)[2] <- "GCFs"
SD_P_F$GCFs[SD_P_F$GCFs == 2] <- "GCF(M)"
SD_P_GCF <- rbind(SD_P_GCF, SD_P_F)
names(SD_C_F)[2] <- "GCFs"
SD_C_F$GCFs[SD_C_F$GCFs == 2] <- "GCF(M)"
SD_C_GCF <- rbind(SD_C_GCF, SD_C_F)
P_N <- merge(N_P_N, SE_P_N, by = c("Target_Type", "NFIX_Target"), all = TRUE) 
C_N <- merge(N_C_N, SE_C_N, by = c("Target_Type", "NFIX_Target"), all = TRUE) 
P_LC <- merge(N_P_LC, SE_P_LC, by = c("Target_Type", "Lifecycle_Target"), all = TRUE) 
C_LC <- merge(N_C_LC, SE_C_LC, by = c("Target_Type", "Lifecycle_Target"), all = TRUE) 
P_FG <- merge(N_P_FG, SE_P_FG, by = c("Target_Type", "FUN_Group_Target"), all = TRUE) 
C_FG <- merge(N_C_FG, SE_C_FG, by = c("Target_Type", "FUN_Group_Target"), all = TRUE) 
P_I <- merge(N_P_I, SE_P_I, by = c("Target_Type", "Type_IND"), all = TRUE) 
C_I <- merge(N_C_I, SE_C_I, by = c("Target_Type", "Type_IND"), all = TRUE) 

#数据输出
##表格
write.xlsx(P_GCF, "result.xlsx", sheetName = "P_GCF")
write.xlsx(C_GCF, "result.xlsx", sheetName = "C_GCF", append = TRUE)
write.xlsx(SD_P_GCF, "result.xlsx", sheetName = "SD_P_GCF", append = TRUE)
write.xlsx(SD_C_GCF, "result.xlsx", sheetName = "SD_C_GCF", append = TRUE)
write.xlsx(SI_P_GCF, "result.xlsx", sheetName = "SI_P_GCF", append = TRUE)
write.xlsx(SI_C_GCF, "result.xlsx", sheetName = "SI_C_GCF", append = TRUE)
write.xlsx(P_N, "result.xlsx", sheetName = "P_N", append = TRUE)
write.xlsx(C_N, "result.xlsx", sheetName = "C_N", append = TRUE)
write.xlsx(P_LC, "result.xlsx", sheetName = "P_LC", append = TRUE)
write.xlsx(C_LC, "result.xlsx", sheetName = "C_LC", append = TRUE)
write.xlsx(P_FG, "result.xlsx", sheetName = "P_FG", append = TRUE)
write.xlsx(C_FG, "result.xlsx", sheetName = "C_FG", append = TRUE)
write.xlsx(P_I, "result.xlsx", sheetName = "P_I", append = TRUE)
write.xlsx(C_I, "result.xlsx", sheetName = "C_I", append = TRUE)


##绘图
###绘图设置
cc <- c("#fb8072", "#ccebc5", "#ccebc5", "#80b1d3", "#fed9a6", "#fed9a6", "#fed9a6", "#fed9a6")
ce <- c("#e41a1c", "#4daf4a", "#4daf4a", "#377eb8", "#ff7f00", "#ff7f00", "#ff7f00", "#ff7f00")
ct <- c("#e41a1c", "#377eb8", "#4daf4a")
c3 <- c("#00afbb", "#e7b800", "#fc4e08")
c7 <- c("#007f87", "#ff5f08", "#00afbb", "#e7b800", "#fc4e08", "#c93e06", "#00d2e1")
#c7 <- c("#1f78b4", "#ff5f08", "#00afbb", "#e7b800", "#fc4e08", "#c93e06", "#6a3d9a")
windowsFonts(myFont = windowsFont("Times New Roman")) #字体修改

###制图
####图1交互关系图
F1a <- ggplot(D1[[1]])+
  geom_col(aes(x = x, y = y, fill = x), width = 0.5)+
  geom_text(aes(x = x, y = y+1.3, label = y), family = "myFont", size = 3)+
  scale_color_manual(values = cc)+
  scale_fill_manual(values = cc)+
  geom_segment(aes(x = 0.75, y = 20, xend = 6.25, yend = 20), color = cc[1], size = 0.7)+
  geom_segment(aes(x = 3.75, y = 37, xend = 6.25, yend = 37), color = cc[4], size = 0.7)+
  annotate("text", label = c("a = 5", "b = 12", "a+b = 17", "a+b = 21", "a+b = -4"), x = c(2, 3, 4, 5, 6), y = c(21.5, 21.5, 21.5, 21.5, 19), colour = ct[1], family = "myFont", fontface = "bold", size = 3)+
  annotate("text", label = c("a*b = 4", "a*b = -16", "AD = CK + a + b"), x = c(5, 6, 4), y = c(38.5, 35.5, 40), colour = ct[2], family = "myFont", fontface = "bold", size = 3)+
  labs(x = "", y = "Response")+
  scale_y_continuous(limits=c(0,44), expand = c(0,0))+
  scale_x_discrete(labels=function(x) str_wrap(x, width=10))+
  theme_bw()+
  theme(legend.position="none",
        title = element_text(family = "myFont",size = 10),
        text = element_text(family = "myFont",size = 7))

F1b <- ggplot(D1[[2]])+
  geom_col(aes(x = x, y = y, fill = x), width = 0.5)+
  geom_text(aes(x = x, y = y+1.3, label = y), family = "myFont", size = 3)+
  scale_color_manual(values = cc)+
  scale_fill_manual(values = cc)+
  geom_segment(aes(x = 0.75, y = 20, xend = 6.25, yend = 20), color = cc[1], size = 0.7)+
  geom_segment(aes(x = 3.75, y = 3, xend = 6.25, yend = 3), color = cc[4], size = 0.7)+
  annotate("text", label = c("a = -5", "b = -12", "a+b = -17", "a+b = -19", "a+b = -8"), x = c(2, 3, 4, 5, 6), y = c(18.5, 18.5, 18.5, 18.5, 18.5), colour = ct[1], family = "myFont", fontface = "bold", size = 3)+
  annotate("text", label = c("a*b = -2", "a*b = 9", "AD = CK + a + b"), x = c(5, 6, 4), y = c(4.5, 4.5, 6), colour = ct[2], family = "myFont", fontface = "bold", size = 3)+
  labs(x = "", y = "Response")+
  scale_y_continuous(limits=c(0,44), expand = c(0,0))+
  scale_x_discrete(labels=function(x) str_wrap(x, width=10))+
  theme_bw()+
  theme(legend.position="none",
        title = element_text(family = "myFont",size = 10),
        text = element_text(family = "myFont",size = 7))

F1c <- ggplot(D1[[3]])+
  geom_col(aes(x = x, y = y, fill = x), width = 0.5)+
  geom_text(aes(x = x, y = y+1.3, label = y), family = "myFont", size = 3)+
  scale_color_manual(values = cc)+
  scale_fill_manual(values = cc)+
  geom_segment(aes(x = 0.75, y = 20, xend = 6.25, yend = 20), color = cc[1], size = 0.7)+
  geom_segment(aes(x = 3.75, y = 13, xend = 6.25, yend = 13), color = cc[4], size = 0.7)+
  annotate("text", label = c("a = 5", "b = -12", "a+b = -17", "a+b = -19", "a+b = -8"), x = c(2, 3, 4, 5, 6), y = c(21.5, 18.5, 18.5, 18.5, 18.5), colour = ct[1], family = "myFont", fontface = "bold", size = 3)+
  annotate("text", label = c("a*b = -2", "a*b = -8", "AD = CK + a + b"), x = c(5, 6, 4), y = c(14.5, 11.5, 16), colour = ct[2], family = "myFont", fontface = "bold", size = 3)+
  labs(x = "", y = "Response")+
  scale_y_continuous(limits=c(0,44), expand = c(0,0))+
  scale_x_discrete(labels=function(x) str_wrap(x, width=10))+
  theme_bw()+
  theme(legend.position="none",
        title = element_text(family = "myFont",size = 10),
        text = element_text(family = "myFont",size = 7))

F1d <- ggplot(D1[[4]])+
  geom_col(aes(x = x, y = y, fill = x), width = 0.4)+
  geom_text(aes(x = x, y = y+1.3, label = y), family = "myFont", size = 3)+
  scale_color_manual(values = cc)+
  scale_fill_manual(values = cc)+
  geom_segment(aes(x = 0.8, y = 20, xend = 4.2, yend = 20), color = cc[1], size = 0.7)+
  geom_segment(aes(x = 3.8, y = 13, xend = 8.2, yend = 13), color = cc[4], size = 0.7)+
  geom_segment(aes(x = 1.8, y = 25, xend = 8.2, yend = 25), color = cc[2], size = 0.7)+
  geom_segment(aes(x = 2.8, y = 8, xend = 8.2, yend = 8), color = cc[3], size = 0.7)+
  annotate("text", label = c("a = 5", "b = -12", "a+b = -7"), x = c(2, 3, 4), y = c(21.5, 18.5, 18.5), colour = ct[1], family = "myFont", fontface = "bold", size = 3)+
  annotate("text", label = c("AD = CK + a + b"), x = c(4), y = c(16), colour = ct[2], family = "myFont", fontface = "bold", size = 3)+
  annotate("text", label = c("a+b > a", "a+b < b", "b < a+b < AD ", "AD < a+b < a "), x = c(5, 6, 7, 8), y = c(27, 6, 9.5, 15), colour = ct[3], family = "myFont", fontface = "bold", size = 3)+
  labs(x = "", y = "Response")+
  scale_y_continuous(limits=c(0,44), expand = c(0,0))+
  scale_x_discrete(labels=function(x) str_wrap(x, width=10))+
  theme_bw()+
  theme(legend.position="none",
        title = element_text(family = "myFont", size = 10),
        text = element_text(family = "myFont", size = 7))

F1e <- ggplot(D1[[5]])+
  geom_segment(aes(x = 4, y = 12, xend = 6, yend = 12), size = 0.5, color = cc[4], linetype = 2)+
  geom_pointrange(aes(x = x, y = y, ymin = y-3, ymax = y+3, color = x))+
  scale_color_manual(values = ce)+
  scale_fill_manual(values = ce)+
  geom_hline(aes(yintercept = 0), linetype = 2)+
  labs(x = "", y = "Response")+
  scale_y_continuous(limits=c(-22,22), expand = c(0,0))+
  scale_x_discrete(labels=function(x) str_wrap(x, width=10))+
  theme_bw()+
  theme(legend.position="none",
        title = element_text(family = "myFont", size = 10),
        text = element_text(family = "myFont", size = 7))

F1f <- ggplot(D1[[6]])+
  geom_segment(aes(x = 4, y = -12, xend = 6, yend = -12), size = 0.5, color = cc[4], linetype = 2)+
  geom_pointrange(aes(x = x, y = y, ymin = y-3, ymax = y+3, color = x))+
  scale_color_manual(values = ce)+
  scale_fill_manual(values = ce)+
  geom_hline(aes(yintercept = 0), linetype = 2)+
  labs(x = "", y = "Response")+
  scale_y_continuous(limits=c(-22,22), expand = c(0,0))+
  scale_x_discrete(labels=function(x) str_wrap(x, width=10))+
  theme_bw()+
  theme(legend.position="none",
        title = element_text(family = "myFont", size = 10),
        text = element_text(family = "myFont", size = 7))

F1g <- ggplot(D1[[7]])+
  geom_rect(aes(xmin = 2, ymin = -7, xmax = 8, ymax = 5), size = 1, alpha = 0.05, fill = "#fff2ae")+
  geom_rect(aes(xmin = 2, ymin = -12, xmax = 8, ymax = -7), size = 1, alpha = 0.05, fill = "#ccebc5")+
  geom_segment(aes(x = 2, y = -7, xend = 8, yend = -7), size = 1, color = cc[4], linetype = 2)+
  geom_pointrange(aes(x = x, y = y, ymin = y-3, ymax = y+3, color = x))+
  scale_color_manual(values = ce)+
  scale_fill_manual(values = ce)+
  geom_hline(aes(yintercept = 0), linetype = 2)+
  labs(x = "", y = "Response")+
  scale_y_continuous(limits=c(-22,22), expand = c(0,0))+
  scale_x_discrete(labels=function(x) str_wrap(x, width=10))+
  theme_bw()+
  theme(legend.position="none",
        title = element_text(family = "myFont", size = 10),
        text = element_text(family = "myFont", size = 7))

####图3数据来源年份与国别
F3a1<- ggplot(D3[[3]][D3[[3]]$Group == "Literature",], aes(x=Year, y=Code, size=Number)) +
  geom_point() +
  theme_bw() +
  theme(legend.position="bottom") +
  scale_y_continuous(breaks=D3[[3]]$Code, labels = D3[[3]]$Country) +
  scale_x_continuous(limits=c(2005,2019), breaks = c(2005,2007,2009,2011,2013,2015,2017,2019)) +
  labs(x = "", y= "") +
  theme(title = element_text(family = "myFont",size=14),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12),
        legend.position ="none",
        panel.grid.minor.y = element_blank(),
        plot.margin=unit(rep(0,4),"lines"))
F3a2 <- ggplot(D3[[2]]) +
  geom_col(aes(x=Year, y=Literature)) +
  scale_x_continuous(expand = c(0.02,0.02)) +
  scale_y_continuous(expand = c(0,0)) +
  theme_void()+
  theme(plot.margin=unit(c(0,0,0,7),"lines"))
F3a3 <- ggplot(D3[[1]]) +
  geom_col(aes(x=Code, y=Literature), width = 0.7) +
  scale_x_continuous(expand = c(0.02,0.02)) +
  scale_y_continuous(expand = c(0,0)) +
  coord_flip() +
  theme_void()+
  theme(plot.margin=unit(c(0,0,1.8,0),"lines"))


F3b1 <- ggplot(D3[[3]][D3[[3]]$Group == "Study",], aes(x=Year, y=Code, size=Number)) +
  geom_point() +
  theme_bw() +
  theme(legend.position="bottom") +
  scale_y_continuous(breaks=D3[[3]]$Code, labels = D3[[3]]$Country) +
  scale_x_continuous(limits=c(2005,2019),breaks = c(2005,2007,2009,2011,2013,2015,2017,2019))+
  labs(x = "", y= "") +
  theme(title=element_text(family = "myFont",size=14),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12),
        legend.position ="none",
        panel.grid.minor.y = element_blank(),
        plot.margin=unit(rep(0,4),"lines"))
F3b2 <- ggplot(D3[[2]]) +
  geom_col(aes(x=Year, y=Study)) +
  scale_x_continuous(expand = c(0.02,0.02)) +
  scale_y_continuous(expand = c(0,0)) +
  theme_void() +
  theme(plot.margin=unit(c(0,0,0,7),"lines"))
F3b3 <- ggplot(D3[[1]]) +
  geom_col(aes(x=Code, y=Study), width = 0.7) +
  scale_x_continuous(expand = c(0.02,0.02)) +
  scale_y_continuous(expand = c(0,0)) +
  coord_flip() +
  theme_void() +
  theme(plot.margin=unit(c(0,0,1.8,0),"lines"))

####图4单因子结果
D4 <- rbind(
  cbind(P_GCF[P_GCF$NUM_GCFs_Record == 1,], 
      data.frame(PC = rep("(A) Performance", length(P_GCF[P_GCF$NUM_GCFs_Record == 1,1])))),
  cbind(C_GCF[C_GCF$NUM_GCFs_Record == 1,], 
        data.frame(PC = rep("(B) Competitiveness", length(C_GCF[C_GCF$NUM_GCFs_Record == 1,1]))))
)
D4$PC <- factor(D4$PC, levels = c("(A) Performance", "(B) Competitiveness"))

F4 <- ggplot(data = D4)+
  geom_point(aes(x = mean, y = GCFs, color = Target_Type, shape = Target_Type),
             size = 3, position = position_dodge(width = -0.5))+
  geom_errorbarh(aes(y = GCFs, color = Target_Type, xmin = down, xmax = up), 
                 height = 0.2, size = 0.8, position = position_dodge(width = -0.5))+
  geom_vline(xintercept = 0, linetype = "dashed", size = 0.7)+
  geom_text(aes(x = up, y = GCFs, color = Target_Type, label = N),
            position=position_dodge(width = -0.5), hjust = -1, show.legend = FALSE)+
  geom_text(data = D4[D4$p < 0.025 | D4$p > 0.975,], aes(x = down, y = GCFs, color = Target_Type), label = "*", size = 8,
            position=position_dodge(width = -0.5), hjust = 1.5, show.legend = FALSE)+
  scale_x_continuous(expand = c(0.1,0.05))+
  scale_color_manual(values = c3)+
  scale_fill_manual(values = c3)+
  theme_bw()+
  scale_y_discrete(limits = factor(c("GCF(S)","C","N","D","P","T"), 
                                   levels=c("GCF(S)","C","N","D","P","T")))+
  labs(x = "Effect size")+
  facet_wrap(vars(PC))+
  theme(axis.title.y=element_blank(),
        legend.background=element_blank(),
        legend.title=element_blank(),
        legend.text=element_text(family = "myFont", size = 12),
        legend.justification=c(0,0), 
        legend.position=c(0.01,0.1),
        title=element_text(family = "myFont",size=14,),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12),
        strip.text.x = element_text(family = "myFont",size=14))

####图5两因子结果
D5 <- rbind(
  cbind(P_GCF[P_GCF$NUM_GCFs_Record == 2,], 
        data.frame(PC = rep("(A) Performance", length(P_GCF[P_GCF$NUM_GCFs_Record == 2,1])))),
  cbind(C_GCF[C_GCF$NUM_GCFs_Record == 2,], 
        data.frame(PC = rep("(B) Competitiveness", length(C_GCF[C_GCF$NUM_GCFs_Record == 2,1]))))
)
D5$PC <- factor(D5$PC, levels = c("(A) Performance", "(B) Competitiveness"))

F5 <- ggplot(data = D5)+
  geom_point(aes(x = mean, y = GCFs, color = Target_Type, shape = Target_Type),
             size = 3, position = position_dodge(width = -0.5))+
  geom_errorbarh(aes(y = GCFs, color = Target_Type, xmin = down, xmax = up), 
                 height = 0.2, size = 0.8, position = position_dodge(width = -0.5))+
  geom_vline(xintercept = 0, linetype = "dashed", size = 0.7)+
  geom_text(aes(x = up, y = GCFs, color = Target_Type, label = N),
            position=position_dodge(width = -0.5), hjust = -1, show.legend = FALSE)+
  geom_text(data = D5[D5$p < 0.025 | D5$p > 0.975,], aes(x = down, y = GCFs, color = Target_Type), label = "*", size = 8,
            position=position_dodge(width = -0.5), hjust = 1.5, show.legend = FALSE)+
  scale_x_continuous(expand = c(0.1,0.05))+
  scale_color_manual(values = c3)+
  scale_fill_manual(values = c3)+
  theme_bw()+
  scale_y_discrete(limits = factor(c("GCF(M)","NC","DC","DN","PN","TC","TN","TD"), 
                                   levels=c("GCF(M)","NC","DC","DN","PN","TC","TN","TD")))+
  labs(x = "Effect size")+
  facet_wrap(vars(PC))+
  theme(axis.title.y=element_blank(),
        legend.background=element_blank(),
        legend.title=element_blank(),
        legend.text=element_text(family = "myFont", size = 12),
        legend.justification=c(0,0), 
        legend.position=c(0.01,0.1),
        title=element_text(family = "myFont",size=14,),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12),
        strip.text.x = element_text(family = "myFont",size=14))

####图6交互结果
D6 <- list(SD_P_GCF, SD_C_GCF, 
           SI_P_GCF[SI_P_GCF$Freq != 0,], 
           SI_C_GCF[SI_C_GCF$Freq != 0,])
D6[[3]] <- rbind(D6[[3]],
                 cbind(GCFs = "GCF(M)",
                 ddply(D6[[3]], c("Target_Type", "inter"), summarize, Freq = round(sum(Freq), 2))))
D6[[3]] <- merge(D6[[3]], 
                 ddply(D6[[3]], c("Target_Type", "GCFs"), summarize, Total = round(sum(Freq), 2)), 
                 by = c("Target_Type", "GCFs"), all = TRUE) 
D6[[3]]$y <- D6[[3]]$Freq/D6[[3]]$Total
D6[[3]]$inter <- factor(D6[[3]]$inter, levels=c("+Synergistic","-Antagonistic","Synergistic","Additive","Antagonistic","+Antagonistic","-Synergistic"))
D6[[3]]$text[D6[[3]]$inter == "+Synergistic"] <- paste("+S: ", round(D6[[3]]$y[D6[[3]]$inter == "+Synergistic"]*100, 2), "%", sep="")
D6[[3]]$text[D6[[3]]$inter == "-Antagonistic"] <- paste("-A: ", round(D6[[3]]$y[D6[[3]]$inter == "-Antagonistic"]*100, 2), "%", sep="")
D6[[3]]$text[D6[[3]]$inter == "Synergistic"] <- paste("S: ", round(D6[[3]]$y[D6[[3]]$inter == "Synergistic"]*100, 2), "%", sep="")
D6[[3]]$text[D6[[3]]$inter == "Additive"] <- paste("AD: ", round(D6[[3]]$y[D6[[3]]$inter == "Additive"]*100, 2), "%", sep="")
D6[[3]]$text[D6[[3]]$inter == "Antagonistic"] <- paste("A: ", round(D6[[3]]$y[D6[[3]]$inter == "Antagonistic"]*100, 2), "%", sep="")
D6[[3]]$text[D6[[3]]$inter == "+Antagonistic"] <- paste("+A: ", round(D6[[3]]$y[D6[[3]]$inter == "+Antagonistic"]*100, 2), "%", sep="")
D6[[3]]$text[D6[[3]]$inter == "-Synergistic"] <- paste("-S: ", round(D6[[3]]$y[D6[[3]]$inter == "-Synergistic"]*100, 2), "%", sep="")
D6[[3]]$text[D6[[3]]$y < 0.05] <- ""
D6[[3]]$x <- paste(D6[[3]]$GCFs, D6[[3]]$Target_Type)
D6[[3]]$x <- factor(D6[[3]]$x, levels = c("GCF(M) Native", "GCF(M) Invasive", "NC Native", "NC Invasive", "DC Invasive", "DC Native", "DN Invasive", "DN Native", "PN Native", "PN Invasive", "TC Native", "TC Invasive", "TN Native", "TN Invasive", "TD Native", "TD Invasive"))

D6[[4]] <- rbind(D6[[4]],
                 cbind(GCFs = "GCF(M)",
                       ddply(D6[[4]], c("Target_Type", "inter"), summarize, Freq = round(sum(Freq), 2))))
D6[[4]] <- merge(D6[[4]], 
                 ddply(D6[[4]], c("Target_Type", "GCFs"), summarize, Total = round(sum(Freq), 2)), 
                 by = c("Target_Type", "GCFs"), all = TRUE) 
D6[[4]]$y <- D6[[4]]$Freq/D6[[4]]$Total
D6[[4]]$inter <- factor(D6[[4]]$inter, levels=c("+Synergistic","-Antagonistic","Synergistic","Additive","Antagonistic","+Antagonistic","-Synergistic"))
D6[[4]]$text[D6[[4]]$inter == "+Synergistic"] <- paste("+S: ", round(D6[[4]]$y[D6[[4]]$inter == "+Synergistic"]*100, 2), "%", sep="")
D6[[4]]$text[D6[[4]]$inter == "-Antagonistic"] <- paste("-A: ", round(D6[[4]]$y[D6[[4]]$inter == "-Antagonistic"]*100, 2), "%", sep="")
D6[[4]]$text[D6[[4]]$inter == "Synergistic"] <- paste("S: ", round(D6[[4]]$y[D6[[4]]$inter == "Synergistic"]*100, 2), "%", sep="")
D6[[4]]$text[D6[[4]]$inter == "Additive"] <- paste("AD: ", round(D6[[4]]$y[D6[[4]]$inter == "Additive"]*100, 2), "%", sep="")
D6[[4]]$text[D6[[4]]$inter == "Antagonistic"] <- paste("A: ", round(D6[[4]]$y[D6[[4]]$inter == "Antagonistic"]*100, 2), "%", sep="")
D6[[4]]$text[D6[[4]]$inter == "+Antagonistic"] <- paste("+A: ", round(D6[[4]]$y[D6[[4]]$inter == "+Antagonistic"]*100, 2), "%", sep="")
D6[[4]]$text[D6[[4]]$inter == "-Synergistic"] <- paste("-S: ", round(D6[[4]]$y[D6[[4]]$inter == "-Synergistic"]*100, 2), "%", sep="")
D6[[4]]$text[D6[[4]]$y < 0.05] <- ""
D6[[4]]$x <- paste(D6[[4]]$GCFs, D6[[4]]$Target_Type)
D6[[4]]$x <- factor(D6[[4]]$x, levels = c("GCF(M) Native", "GCF(M) Invasive", "NC Native", "NC Invasive", "DC Invasive", "DC Native", "DN Invasive", "DN Native", "PN Native", "PN Invasive", "TC Native", "TC Invasive", "TN Native", "TN Invasive", "TD Native", "TD Invasive"))


F6a1 <- ggplot(D6[[1]])+
  geom_point(aes(x = mean, y = GCFs, color = Target_Type, shape = Target_Type),
             size = 3, position = position_dodge(width = -0.5))+
  geom_errorbarh(aes(y = GCFs, color = Target_Type, xmin = down, xmax = up), 
                 height = 0.2, size = 0.8, position = position_dodge(width = -0.5))+
  geom_vline(xintercept = 0, linetype = "dashed", size = 0.7)+
  geom_text(data = D6[[1]][D6[[1]]$p < 0.025 | D6[[1]]$p > 0.975,], aes(x = down, y = GCFs, color = Target_Type), label = "*", size = 8,
            position=position_dodge(width = -0.5), hjust = 1.5, show.legend = FALSE)+
  #scale_x_continuous(expand = c(0.1,0.05))+
  scale_color_manual(values = c3)+
  scale_fill_manual(values = c3)+
  theme_bw()+
  scale_y_discrete(limits = factor(c("GCF(M)","NC","DC","DN","PN","TC","TN","TD"), 
                                   levels=c("GCF(M)","NC","DC","DN","PN","TC","TN","TD")))+
  labs(x = "Difference size")+
  theme(axis.title.y=element_blank(),
        legend.background=element_blank(),
        legend.title=element_blank(),
        legend.text=element_text(family = "myFont", size = 14),
        title=element_text(family = "myFont",size=14,),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12),
        strip.text.x = element_text(family = "myFont",size=14))

F6b1 <- ggplot(D6[[2]])+
  geom_point(aes(x = mean, y = GCFs, color = Target_Type, shape = Target_Type),
             size = 3, position = position_dodge(width = -0.5))+
  geom_errorbarh(aes(y = GCFs, color = Target_Type, xmin = down, xmax = up), 
                 height = 0.2, size = 0.8, position = position_dodge(width = -0.5))+
  geom_vline(xintercept = 0, linetype = "dashed", size = 0.7)+
  geom_text(data = D6[[2]][D6[[2]]$p < 0.025 | D6[[2]]$p > 0.975,], aes(x = down, y = GCFs, color = Target_Type), label = "*", size = 8,
            position=position_dodge(width = -0.5), hjust = 1.5, show.legend = FALSE)+
  scale_x_continuous(expand = c(0.1,0.05))+
  scale_color_manual(values = c3)+
  scale_fill_manual(values = c3)+
  theme_bw()+
  scale_y_discrete(limits = factor(c("GCF(M)","NC","DC","DN","PN","TC","TN","TD"), 
                                   levels=c("GCF(M)","NC","DC","DN","PN","TC","TN","TD")))+
  labs(x = "Difference size")+
  theme(axis.title.y=element_blank(),
        legend.background=element_blank(),
        legend.title=element_blank(),
        legend.text=element_text(family = "myFont", size = 14),
        title=element_text(family = "myFont",size=14,),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12),
        strip.text.x = element_text(family = "myFont",size=14))

F6a2 <- ggplot(D6[[3]])+
  geom_col(aes(x = x, y = y, fill = inter, color = inter), width = 0.6)+
  geom_text(aes(x = x, y = y, label = text, fill = inter), position = position_stack(vjust = 0.5), color = "#FFFFFF", size = 3, fontface = "bold")+
  scale_y_continuous(expand = c(0,0))+
  scale_color_manual(values = c7)+
  scale_fill_manual(values = c7)+
  coord_flip()+
  theme_void()+
  theme(plot.margin=unit(c(0,0,2,0),"lines"),
        legend.title=element_blank(),
        legend.text=element_text(family = "myFont", size = 8),
        legend.key.size =  unit(1, "lines"))

F6b2 <- ggplot(D6[[4]])+
  geom_col(aes(x = x, y = y, fill = inter, color = inter), width = 0.6)+
  geom_text(aes(x = x, y = y, label = text, fill = inter), position = position_stack(vjust = 0.5), color = "#FFFFFF", size = 3, fontface = "bold", show.legend = FALSE)+
  scale_y_continuous(expand = c(0,0))+
  scale_color_manual(values = c7)+
  scale_fill_manual(values = c7)+
  coord_flip()+
  theme_void()+
  theme(plot.margin=unit(c(0,0,2,0),"lines"),
        legend.title=element_blank(),
        legend.text=element_text(family = "myFont", size = 8),
        legend.key.size =  unit(1, "lines"))


####图S2其它分组图
DS2 <- list(rbind(cbind(P_N, data.frame(PC = rep("(A) Performance", length(P_N[,1])))),
                  cbind(C_N, data.frame(PC = rep("(B) Competitiveness", length(C_N[,1]))))),
            rbind(cbind(P_LC, data.frame(PC = rep("(A) Performance", length(P_LC[,1])))),
                  cbind(C_LC, data.frame(PC = rep("(B) Competitiveness", length(C_LC[,1]))))),
            rbind(cbind(P_FG, data.frame(PC = rep("(A) Performance", length(P_FG[,1])))),
                  cbind(C_FG, data.frame(PC = rep("(B) Competitiveness", length(C_FG[,1]))))),
            rbind(cbind(P_I, data.frame(PC = rep("(A) Performance", length(P_I[,1])))),
                  cbind(C_I, data.frame(PC = rep("(B) Competitiveness", length(C_I[,1]))))))
DS2.1 <- c(names(DS2[[1]])[2], names(DS2[[2]])[2], names(DS2[[3]])[2], names(DS2[[4]])[2])
FS2 <- list()

for (i in c(1:4)) {
  FS2[[i]] <- ggplot(DS2[[i]])+
    #geom_point(aes_string(x = "mean", y = DS2.1[i], color = "Target_Type", shape = "Target_Type"),
    #           size = 3)+
    #geom_errorbarh(aes_string(y = DS2.1[i], color = "Target_Type", xmin = "down", xmax = "up"), 
    #height = 0.2, size = 0.8)+
    geom_pointrange(aes_string(x = "mean", y = DS2.1[i], color = "Target_Type", xmin = "down", xmax = "up"), 
                    fatten = 2, size = 1, position=position_dodge(width = -0.5))+
    geom_vline(xintercept = 0, linetype = "dashed", size = 0.7)+
    geom_text(aes_string(x = "up", y = DS2.1[i], color = "Target_Type", label = "N"),
               hjust = -0.5, show.legend = FALSE, position=position_dodge(width = -0.5))+
    #geom_text(data = DS2[[i]][DS2[[i]]$p < 0.025 | DS2[[i]]$p > 0.975,], 
    #          aes_string(x = "down", y = DS2.1[i], color =" Target_Type"), 
    #          label = "*", size = 8, position=position_dodge(width = -0.5),
    #          hjust = 1.5, show.legend = FALSE)+
    scale_x_continuous(expand = c(0.2,0.2))+
    scale_color_manual(values = c3)+
    scale_fill_manual(values = c3)+
    theme_bw()+
    labs(x = "Effect size")+
    facet_wrap(vars(PC), scales= "free_x")+
    theme(axis.title.y=element_blank(),
          legend.background=element_blank(),
          legend.title=element_blank(),
          legend.text=element_text(family = "myFont", size = 12),
          legend.justification=c(0,0), 
          legend.position=c(0.01,0.01),
          title=element_text(family = "myFont",size=14),
          axis.text.x = element_text(family = "myFont",size=12),
          axis.text.y = element_text(family = "myFont",size=12),
          strip.text.x = element_text(family = "myFont",size=14))
}


###合图
F1 <- ggarrange(
  ggarrange(F1a,F1b,F1c, nrow = 3, ncol = 1, labels = LETTERS[1:3], font.label = list(color = "black")),
  ggarrange(F1d, 
            ggarrange(F1e,F1f,nrow = 1, ncol = 2, labels = LETTERS[5:6], font.label = list(color = "black")),
            F1g, nrow = 3, ncol = 1, labels = c(LETTERS[4], "", LETTERS[7]), font.label = list(color = "black")),
  nrow = 1, ncol = 2, widths = c(1,1.6)
)

F3 <- ggarrange(
  ggarrange(F3a2, NULL, F3a1, F3a3, nrow = 2, ncol = 2, widths = c(8,1), heights = c(1,8),common.legend = T, legend = "left"),
  ggarrange(F3b2, NULL, F3b1, F3b3, nrow = 2, ncol = 2, widths = c(8,1), heights = c(1,8),common.legend = T, legend = "left"),
  nrow = 2, ncol = 1, labels = c(LETTERS[1:2])
)

F6 <- ggarrange(
  ggarrange(F6a1, F6b1, nrow = 2, ncol = 1, common.legend = T, legend = "top", labels = c(LETTERS[1:2])),
  ggarrange(F6a2, F6b2, nrow = 2, ncol = 1, common.legend = T, legend = "top"),
  nrow = 1, ncol = 2, align = "v", widths = c(1.5,1))

FS2T <- ggarrange(
  ggarrange(NULL, FS2[[1]], NULL, FS2[[3]], nrow = 4, ncol = 1, heights = c(0.3,2,0.2,6),
            hjust = 0,  vjust = -0.2, align = "v", common.legend = T, legend = "top",
            labels = c("", "(A) Nitrogen fixation", "", "(C) Functional group")),
  ggarrange(NULL, FS2[[2]], NULL, FS2[[4]], nrow = 4, ncol = 1, heights = c(0.3,4,0.2,3),
            hjust = 0,  vjust = -0.2, align = "v", common.legend = T, legend = "none",
            labels = c("", "(B) Life cycle", "", "(D) Performance response")),
  nrow = 1, ncol = 2, common.legend = T)


###存图
ggsave(file = "F1.svg",plot = F1, width = 10, height = 10)
ggsave(file = "F3.svg",plot = F3, width = 10, height = 8)
ggsave(file = "F4.svg",plot = F4, width = 8, height = 6)
ggsave(file = "F5.svg",plot = F5, width = 8, height = 7.5)
ggsave(file = "F6.svg",plot = F6, width = 10, height = 10)
ggsave(file = "FS2.svg",plot = FS2T, width = 12, height = 10)

